##################################################################################
######## Replication: Border Barriers and Local Armed Uprisings ##################
######## Security Studies, 2025 ##################################################
######## Md Muhibbur Rahman ######################################################
######## Cornell University ######################################################
######## March, 2025 #############################################################
##################################################################################

########################## Causal mediation analysis

### Load Libraries

library(ggplot2)
library(broom)
library(dplyr)
library(mediation)
library(haven)

### Load data


# Load data
data <- read.dta13("ss_replication_data.dta")



###################################################################################
###### Figure 5:Causal mediation effects (GDP per capita and mountainous terrain)
##################################################################################


### Running the Models

# State capacity measure 1: combination of GDP per capita and rugged terrain

# Ensure build_cap_main is treated as a factor
data$build_cap_main <- as.factor(data$build_cap_main)
# build_cap_main a categorical variable with three categories, 1: low GDP pc and high rugged terrain
# GDP pc high and low is defined whether country is more than mean+1SD
# Rugged terrain high and low is defined whether country is more than mean+1SD 



# Mediator model (binary logistic regression)
mediator_model <- glm(bw_dummy ~ build_cap_main + riv_dummy + threat + neighbors_number + prev_civ_conf + extek + oilpcl + polity2_P4 + lnpop + postcold +
                        region_SSA + region_SoAsia + region_MENA + region_Latin + region_EAsiaPac + region_EurCentAsia + region_NoAmer, 
                      data = data, family = binomial)

summary(mediator_model)


# Outcome model (binary logistic regression)
outcome_model <- glm(civ_war ~ build_cap_main + bw_dummy + riv_dummy + threat + neighbors_number + prev_civ_conf + extek + oilpcl + polity2_P4 + lnpop + postcold +
                       region_SSA + region_SoAsia + region_MENA + region_Latin + region_EAsiaPac + region_EurCentAsia + region_NoAmer, 
                     data = data, family = binomial)

summary(outcome_model)

# Perform mediation analysis
mediation_result <- mediate(mediator_model, outcome_model, treat = "build_cap_main", mediator = "bw_dummy", boot = TRUE, sims = 1000)

# Summarize the mediation analysis results
summary(mediation_result)



### Figure 5 (Upper Panel)

# Tidy the mediator model coefficients
mediator_coeffs <- tidy(mediator_model)

# Tidy the outcome model coefficients
outcome_coeffs <- tidy(outcome_model)

# Combine the two datasets
mediator_coeffs$model <- "Mediator Model"
outcome_coeffs$model <- "Outcome Model"

coefficients <- bind_rows(mediator_coeffs, outcome_coeffs)

# Filter out the unwanted region variables and "mtnest"
coefficients_filtered <- subset(coefficients, 
                                !term %in% c("region_SSA", "region_SoAsia", "region_MENA", 
                                             "region_Latin", "region_EAsiaPac", 
                                             "region_EurCentAsia", "region_NoAmer", "mtnest"))

# Rename the terms, including renaming the intercept
coefficients_filtered$term <- recode(coefficients_filtered$term,
                                     "(Intercept)" = "Constant",
                                     "bw_dummy" = "Border barrier",
                                     "riv_dummy" = "Interstate rivalry",
                                     "threat" = "External threat",
                                     "neighbors_number" = "Number of neighbors",
                                     "prev_civ_conf" = "Number of previous civil wars",
                                     "extek" = "Excluded group with large TEK",
                                     "oilpcl" = "Oil production per capita",
                                     "polity2_P4" = "Democracy",
                                     "lnpop" = "Total population (logged)",
                                     "build_cap_main2" = "Medium state capacity",
                                     "build_cap_main3" = "High state capacity",
                                     "postcold" = "Cold War"
)

# Order the terms in the specified order (without Mountainous terrain)
order_of_terms <- c(
  "Medium state capacity", "High state capacity", "Border barrier", 
  "Interstate rivalry", "External threat", "Number of neighbors", 
  "Number of previous civil wars", "Excluded group with large TEK", 
  "Oil production per capita", "Democracy", "Total population (logged)", 
  "Cold War", "Constant"
)
coefficients_filtered$term <- factor(coefficients_filtered$term, levels = rev(order_of_terms))

# Plot the filtered, renamed, and ordered coefficients with a reference line at x = 0
ggplot(coefficients_filtered, aes(x = estimate, y = term, color = model)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = estimate - std.error, xmax = estimate + std.error), height = 0.2) +
  facet_wrap(~model, scales = "free") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +  # Add reference line at x = 0
  labs(title = "",
       x = "Estimate",
       y = "Term") +
  theme_minimal()

# Save as TIFF file 
ggsave("Coefficient_Plot.tiff", plot = last_plot(), dpi = 300, width = 8, height = 4, units = "in", device = "tiff", bg = "white")



### Figure 5 (Bottom Panel)


# Plot the ACME and ADE for the average effects with custom y-axis titles
acme_ade_average2 <- data.frame(
  effect = c("Average Causal Mediation Effect", "Average Direct Effect"),
  estimate = c(mediation_result$d.avg, mediation_result$z.avg),
  lower = c(mediation_result$d.avg.ci[1], mediation_result$z.avg.ci[1]),
  upper = c(mediation_result$d.avg.ci[2], mediation_result$z.avg.ci[2])
)

# Create plot
acme_ade_plot <- ggplot(acme_ade_average2, aes(x = estimate, y = effect)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +  # Reference line at x = 0
  labs(x = "Estimate", y = "") +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "white", color = NA), # Ensure white background
    plot.background = element_rect(fill = "white", color = NA)
  )

# Save as  TIFF file 
ggsave("ACME_ADE_Plot.tiff", plot = acme_ade_plot, dpi = 300, width = 8, height = 2, units = "in", device = "tiff", bg = "white")





###################################################################################
###### Figure A5: Causal mediation effects (H&S measure and mountainous terrain) 
##################################################################################




### Running the Models

# State capacity measure 2: combination of Hanson and Sigman's state capacity measure and rugged terrain


# Ensure build_cap_main is treated as a factor
data$build_cap_alter <- as.factor(data$build_cap_alter)
# build_cap_main2 a categorical variable with three categories, 1: low GDP pc and high rugged terrain
# GDP pc high and low is defined whether country is more than mean+1SD
# Rugged terrain high and low is defined whether country is more than mean+1SD 



# Mediator model (binary logistic regression)
mediator_model2 <- glm(bw_dummy ~ build_cap_alter + riv_dummy + threat + neighbors_number + prev_civ_conf + extek + oilpcl + polity2_P4 + lnpop + postcold +
                         region_SSA + region_SoAsia + region_MENA + region_Latin + region_EAsiaPac + region_EurCentAsia + region_NoAmer, 
                       data = data, family = binomial)

summary(mediator_model2)

# Outcome model (binary logistic regression)
outcome_model2 <- glm(civ_war ~ build_cap_alter + bw_dummy + riv_dummy + threat + neighbors_number + prev_civ_conf + extek + oilpcl + polity2_P4 + lnpop + postcold +
                        region_SSA + region_SoAsia + region_MENA + region_Latin + region_EAsiaPac + region_EurCentAsia + region_NoAmer, 
                      data = data, family = binomial)

summary(outcome_model2)

# Perform mediation analysis
mediation_result2 <- mediate(mediator_model2, outcome_model2, treat = "build_cap_alter", mediator = "bw_dummy", boot = TRUE, sims = 1000)

# Summarize the mediation analysis results
summary(mediation_result2)




### Figure A5 (Upper Panel)


# Tidy the mediator model coefficients
mediator_coeffs2 <- tidy(mediator_model2)

# Tidy the outcome model coefficients
outcome_coeffs2 <- tidy(outcome_model2)

# Combine the two datasets
mediator_coeffs2$model <- "Mediator Model2"
outcome_coeffs2$model <- "Outcome Model2"

coefficients2 <- bind_rows(mediator_coeffs2, outcome_coeffs2)

# Filter out the unwanted region variables and "mtnest"
coefficients_filtered2 <- subset(coefficients2, 
                                 !term %in% c("region_SSA", "region_SoAsia", "region_MENA", 
                                              "region_Latin", "region_EAsiaPac", 
                                              "region_EurCentAsia", "region_NoAmer", "mtnest"))

# Rename the terms, including renaming the intercept
coefficients_filtered2$term <- recode(coefficients_filtered2$term,
                                      "(Intercept)" = "Constant",
                                      "bw_dummy" = "Border barrier",
                                      "riv_dummy" = "Interstate rivalry",
                                      "threat" = "External threat",
                                      "neighbors_number" = "Number of neighbors",
                                      "prev_civ_conf" = "Number of previous civil wars",
                                      "extek" = "Excluded group with large TEK",
                                      "oilpcl" = "Oil production per capita",
                                      "polity2_P4" = "Democracy",
                                      "lnpop" = "Total population (logged)",
                                      "build_cap_alter2" = "Medium state capacity",
                                      "build_cap_alter3" = "High state capacity",
                                      "postcold" = "Cold War"
)

# Order the terms in the specified order (without Mountainous terrain)
order_of_terms2 <- c(
  "Medium state capacity", "High state capacity", "Border barrier", 
  "Interstate rivalry", "External threat", "Number of neighbors", 
  "Number of previous civil wars", "Excluded group with large TEK", 
  "Oil production per capita", "Democracy", "Total population (logged)", 
  "Cold War", "Constant"
)
coefficients_filtered2$term <- factor(coefficients_filtered2$term, levels = rev(order_of_terms2))

# Plot the filtered, renamed, and ordered coefficients with a reference line at x = 0
ggplot(coefficients_filtered2, aes(x = estimate, y = term, color = model)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = estimate - std.error, xmax = estimate + std.error), height = 0.2) +
  facet_wrap(~model, scales = "free") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +  # Add reference line at x = 0
  labs(title = "",
       x = "Estimate",
       y = "Term") +
  theme_minimal()


# Save as high-resolution TIFF file (300 DPI, 2400x1200 pixels)
ggsave("Coefficient_Plot2.tiff", plot = last_plot(), dpi = 300, width = 8, height = 4, units = "in", device = "tiff", bg = "white")




### Figure A5 (Bottom Panel)


## Plot ADE AND ACME
acme_ade_average2 <- data.frame(
  effect = c("Average Causal Mediation Effect", "Average Direct Effect"),
  estimate = c(mediation_result2$d.avg, mediation_result2$z.avg),
  lower = c(mediation_result2$d.avg.ci[1], mediation_result2$z.avg.ci[1]),
  upper = c(mediation_result2$d.avg.ci[2], mediation_result2$z.avg.ci[2])
)


acme_ade_plot2 <- ggplot(acme_ade_average2, aes(x = estimate, y = effect)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +  # Add reference line at x = 0
  labs(title = "",
       x = "Estimate",
       y = "") +
  theme_minimal()+
  theme(
    panel.background = element_rect(fill = "white", color = NA), # Ensure white background
    plot.background = element_rect(fill = "white", color = NA)
  )

# Save as TIFF file 
ggsave("ACME_ADE_Plot2.tiff", plot = acme_ade_plot2, dpi = 300, width = 8, height = 2, units = "in", device = "tiff", bg = "white")






###########################################################################################
#################### Robustness checks with extended data #################################
###########################################################################################



data <- read.dta13("ss_replication_data.dta")



###################################################################################################
###### Figure A6:  Causal mediation effects (GDP per capita and mountainous terrain, extended data) 
###################################################################################################



### Running the Models 

# State capacity measure 1: combination of GDP per capita and rugged terrain

# Ensure build_cap_main is treated as a factor
data$build_cap_main <- as.factor(data$build_cap_main)
# build_cap_main a categorical variable with three categories, 1: low GDP pc and high rugged terrain
# GDP pc high and low is defined whether country is more than mean+1SD
# Rugged terrain high and low is defined whether country is more than mean+1SD 


# Mediator model (binary logistic regression)
mediator_model3 <- glm(bw_extended_KPS ~ build_cap_main + riv_dummy + threat + neighbors_number + prev_civ_conf + extek + oilpcl + polity2_P4 + lnpop, 
                       data = data, family = binomial)

summary(mediator_model3)


# Outcome model (binary logistic regression)
outcome_model3 <- glm(civ_war ~ build_cap_main + bw_extended_KPS + riv_dummy + threat + neighbors_number + prev_civ_conf + extek + oilpcl + polity2_P4 + lnpop, 
                      data = data, family = binomial)

summary(outcome_model3)

# Perform mediation analysis
mediation_result3 <- mediate(mediator_model3, outcome_model3, treat = "build_cap_main", mediator = "bw_extended_KPS", boot = TRUE, sims = 1000)

# Summarize the mediation analysis results
summary(mediation_result3)



### Figure A6 (Upper Panel)

# Tidy the mediator model coefficients
mediator_coeffs3 <- tidy(mediator_model3)

# Tidy the outcome model coefficients
outcome_coeffs3 <- tidy(outcome_model3)

# Combine the two datasets
mediator_coeffs3$model <- "Mediator Model3"
outcome_coeffs3$model <- "Outcome Model3"

coefficients3 <- bind_rows(mediator_coeffs3, outcome_coeffs3)


# Rename the terms, including renaming the intercept
coefficients3$term <- recode(coefficients3$term,
                             "(Intercept)" = "Constant",
                             "bw_extended_KPS" = "Border barrier (extended)",
                             "riv_dummy" = "Interstate rivalry",
                             "threat" = "External threat",
                             "neighbors_number" = "Number of neighbors",
                             "prev_civ_conf" = "Number of previous civil wars",
                             "extek" = "Excluded group with large TEK",
                             "oilpcl" = "Oil production per capita",
                             "polity2_P4" = "Democracy",
                             "lnpop" = "Total population (logged)",
                             "build_cap_main2" = "Medium state capacity",
                             "build_cap_main3" = "High state capacity"
)


# Order the terms in the specified order (without Mountainous terrain)
order_of_terms3 <- c(
  "Medium state capacity", "High state capacity", "Border barrier (extended)", 
  "Interstate rivalry", "External threat", "Number of neighbors", 
  "Number of previous civil wars", "Excluded group with large TEK", 
  "Oil production per capita", "Democracy", "Total population (logged)", "Constant"
)
coefficients3$term <- factor(coefficients3$term, levels = rev(order_of_terms3))

# Plot the filtered, renamed, and ordered coefficients with a reference line at x = 0
ggplot(coefficients3, aes(x = estimate, y = term, color = model)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = estimate - std.error, xmax = estimate + std.error), height = 0.2) +
  facet_wrap(~model, scales = "free") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +  # Add reference line at x = 0
  labs(title = "",
       x = "Estimate",
       y = "Term") +
  theme_minimal()



# Save as high-resolution TIFF file (300 DPI, 2400x1200 pixels)
ggsave("Coefficient_Plot3.tiff", plot = last_plot(), dpi = 300, width = 8, height = 4, units = "in", device = "tiff", bg = "white")



### Figure A6 (Bottom Panel)


# Plot the ACME and ADE for the average effects with custom y-axis titles

acme_ade_average3 <- data.frame(
  effect = c("Average Causal Mediation Effect", "Average Direct Effect"),
  estimate = c(mediation_result3$d.avg, mediation_result3$z.avg),
  lower = c(mediation_result3$d.avg.ci[1], mediation_result3$z.avg.ci[1]),
  upper = c(mediation_result3$d.avg.ci[2], mediation_result3$z.avg.ci[2])
)

acme_ade_plot3 <- ggplot(acme_ade_average3, aes(x = estimate, y = effect)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +  # Add reference line at x = 0
  labs(title = "",
       x = "Estimate",
       y = "") +
  theme_minimal()+
  theme(
    panel.background = element_rect(fill = "white", color = NA), # Ensure white background
    plot.background = element_rect(fill = "white", color = NA)
  )

# Save as high-resolution TIFF file (300 DPI)
ggsave("ACME_ADE_Plot3.tiff", plot = acme_ade_plot3, dpi = 300, width = 8, height = 2, units = "in", device = "tiff", bg = "white")





#####################################################################################################
###### Figure A7:  Causal mediation effects (the H&S measure and mountainous terrain, extended data) 
#####################################################################################################


### Running the Models

# State capacity measure 1: combination of GDP per capita and rugged terrain

# Ensure build_cap_main is treated as a factor
data$build_cap_alter <- as.factor(data$build_cap_alter)
# build_cap_main a categorical variable with three categories, 1: low GDP pc and high rugged terrain
# GDP pc high and low is defined whether country is more than mean+1SD
# Rugged terrain high and low is defined whether country is more than mean+1SD 


# Mediator model (binary logistic regression)
mediator_model4 <- glm(bw_extended_KPS ~ build_cap_alter + riv_dummy + threat + neighbors_number + prev_civ_conf + extek + oilpcl + polity2_P4 + lnpop, 
                       data = data, family = binomial)

summary(mediator_model4)


# Outcome model (binary logistic regression)
outcome_model4 <- glm(civ_war ~ build_cap_alter + bw_extended_KPS + riv_dummy + threat + neighbors_number + prev_civ_conf + extek + oilpcl + polity2_P4 + lnpop, 
                      data = data, family = binomial)

summary(outcome_model4)

# Perform mediation analysis
mediation_result4 <- mediate(mediator_model4, outcome_model4, treat = "build_cap_alter", mediator = "bw_extended_KPS", boot = TRUE, sims = 1000)

# Summarize the mediation analysis results
summary(mediation_result4)



### Figure 7 (Upper Panel)

# Tidy the mediator model coefficients
mediator_coeffs4 <- tidy(mediator_model4)

# Tidy the outcome model coefficients
outcome_coeffs4 <- tidy(outcome_model4)

# Combine the two datasets
mediator_coeffs4$model <- "Mediator Model4"
outcome_coeffs4$model <- "Outcome Model4"

coefficients4 <- bind_rows(mediator_coeffs4, outcome_coeffs4)


# Rename the terms, including renaming the intercept
coefficients4$term <- recode(coefficients4$term,
                             "(Intercept)" = "Constant",
                             "bw_extended_KPS" = "Border barrier (extended)",
                             "riv_dummy" = "Interstate rivalry",
                             "threat" = "External threat",
                             "neighbors_number" = "Number of neighbors",
                             "prev_civ_conf" = "Number of previous civil wars",
                             "extek" = "Excluded group with large TEK",
                             "oilpcl" = "Oil production per capita",
                             "polity2_P4" = "Democracy",
                             "lnpop" = "Total population (logged)",
                             "build_cap_alter2" = "Medium state capacity",
                             "build_cap_alter3" = "High state capacity"
)


# Order the terms in the specified order (without Mountainous terrain)
order_of_terms4 <- c(
  "Medium state capacity", "High state capacity", "Border barrier (extended)", 
  "Interstate rivalry", "External threat", "Number of neighbors", 
  "Number of previous civil wars", "Excluded group with large TEK", 
  "Oil production per capita", "Democracy", "Total population (logged)", "Constant"
)
coefficients4$term <- factor(coefficients4$term, levels = rev(order_of_terms4))

# Plot the filtered, renamed, and ordered coefficients with a reference line at x = 0
ggplot(coefficients4, aes(x = estimate, y = term, color = model)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = estimate - std.error, xmax = estimate + std.error), height = 0.2) +
  facet_wrap(~model, scales = "free") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +  # Add reference line at x = 0
  labs(title = "",
       x = "Estimate",
       y = "Term") +
  theme_minimal()


# Save as high-resolution TIFF file (300 DPI, 2400x1200 pixels)
ggsave("Coefficient_Plot4.tiff", plot = last_plot(), dpi = 300, width = 8, height = 4, units = "in", device = "tiff", bg = "white")


### Figure 7 (Bottom Panel)


# Plot the ACME and ADE for the average effects with custom y-axis titles

acme_ade_average4 <- data.frame(
  effect = c("Average Causal Mediation Effect", "Average Direct Effect"),
  estimate = c(mediation_result4$d.avg, mediation_result4$z.avg),
  lower = c(mediation_result4$d.avg.ci[1], mediation_result4$z.avg.ci[1]),
  upper = c(mediation_result4$d.avg.ci[2], mediation_result4$z.avg.ci[2])
)


acme_ade_plot4 <- ggplot(acme_ade_average4, aes(x = estimate, y = effect)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +  # Add reference line at x = 0
  labs(title = "",
       x = "Estimate",
       y = "") +
  theme_minimal()+
  theme(
    panel.background = element_rect(fill = "white", color = NA), # Ensure white background
    plot.background = element_rect(fill = "white", color = NA)
  )

# Save as high-resolution TIFF file (300 DPI)
ggsave("ACME_ADE_Plot4.tiff", plot = acme_ade_plot4, dpi = 300, width = 8, height = 2, units = "in", device = "tiff", bg = "white")






###################################################################################
#################### Additional Robustness checks #################################
###################################################################################


data <- read.dta13("ss_replication_data.dta")



##################################################################
###### Figure A8:  Causal mediation effects (GDP per capita only) 
##################################################################



## State capacity measure 2: GDP per capita only


# Mediator model (binary logistic regression)
mediator_model5 <- glm(bw_dummy ~ lngdppc + riv_dummy + threat + neighbors_number + prev_civ_conf + extek + mtnest + oilpcl + polity2_P4 + lnpop + postcold +
                         region_SSA + region_SoAsia + region_MENA + region_Latin + region_EAsiaPac + region_EurCentAsia + region_NoAmer, 
                       data = data, family = binomial)

summary(mediator_model5)



# Outcome model (binary logistic regression)
outcome_model5 <- glm(civ_war ~ lngdppc + bw_dummy + riv_dummy + threat + neighbors_number + prev_civ_conf + extek + +mtnest + oilpcl + polity2_P4 + lnpop + postcold +
                        region_SSA + region_SoAsia + region_MENA + region_Latin + region_EAsiaPac + region_EurCentAsia + region_NoAmer, 
                      data = data, family = binomial)

summary(outcome_model5)

# Perform mediation analysis
mediation_result <- mediate(mediator_model5, outcome_model5, treat = "lngdppc", mediator = "bw_dummy", boot = TRUE, sims = 1000)

# Summarize the mediation analysis results
summary(mediation_result)



### Figure 8 (Upper Panel)


# Tidy the mediator model coefficients
mediator_coeffs <- tidy(mediator_model5)

# Tidy the outcome model coefficients
outcome_coeffs <- tidy(outcome_model5)

# Combine the two datasets
mediator_coeffs$model <- "Mediator Model5"
outcome_coeffs$model <- "Outcome Model5"

coefficients <- bind_rows(mediator_coeffs, outcome_coeffs)

# Filter out the unwanted region variables and "mtnest"
coefficients_filtered <- subset(coefficients, 
                                !term %in% c("region_SSA", "region_SoAsia", "region_MENA", 
                                             "region_Latin", "region_EAsiaPac", 
                                             "region_EurCentAsia", "region_NoAmer", "mtnest"))

# Rename the terms, including renaming the intercept
coefficients_filtered$term <- recode(coefficients_filtered$term,
                                     "(Intercept)" = "Constant",
                                     "bw_dummy" = "Border barrier",
                                     "riv_dummy" = "Interstate rivalry",
                                     "threat" = "External threat",
                                     "neighbors_number" = "Number of neighbors",
                                     "prev_civ_conf" = "Number of previous civil wars",
                                     "extek" = "Excluded group with large TEK",
                                     "mtnest"= "Mountainous terrain",
                                     "oilpcl" = "Oil production per capita",
                                     "polity2_P4" = "Democracy",
                                     "lnpop" = "Total population (logged)",
                                     "lngdppc" = "GDP per capita (logged)",
                                     "postcold" = "Cold War"
)

# Order the terms in the specified order (without Mountainous terrain)
order_of_terms <- c(
  "GDP per capita (logged)", "Border barrier", "Interstate rivalry", "External threat", "Number of neighbors", 
  "Number of previous civil wars", "Excluded group with large TEK", "Mountainous terrain",
  "Oil production per capita", "Democracy", "Total population (logged)", 
  "Cold War", "Constant"
)
coefficients_filtered$term <- factor(coefficients_filtered$term, levels = rev(order_of_terms))

# Plot the filtered, renamed, and ordered coefficients with a reference line at x = 0
ggplot(coefficients_filtered, aes(x = estimate, y = term, color = model)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = estimate - std.error, xmax = estimate + std.error), height = 0.2) +
  facet_wrap(~model, scales = "free") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +  # Add reference line at x = 0
  labs(title = "",
       x = "Estimate",
       y = "Term") +
  theme_minimal()



# Save as high-resolution TIFF file (300 DPI, 2400x1200 pixels)
ggsave("Coefficient_Plot5.tiff", plot = last_plot(), dpi = 300, width = 8, height = 4, units = "in", device = "tiff", bg = "white")



### Figure 8 (Bottom Panel)

# Plot the ACME and ADE for the average effects with custom y-axis titles

acme_ade_average5 <- data.frame(
  effect = c("Average Causal Mediation Effect", "Average Direct Effect"),
  estimate = c(mediation_result$d.avg, mediation_result$z.avg),
  lower = c(mediation_result$d.avg.ci[1], mediation_result$z.avg.ci[1]),
  upper = c(mediation_result$d.avg.ci[2], mediation_result$z.avg.ci[2])
)


acme_ade_plot5 <- ggplot(acme_ade_average5, aes(x = estimate, y = effect)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +  # Add reference line at x = 0
  labs(title = "",
       x = "Estimate",
       y = "") +
  theme_minimal()+
  theme(
    panel.background = element_rect(fill = "white", color = NA), # Ensure white background
    plot.background = element_rect(fill = "white", color = NA)
  )

# Save as high-resolution TIFF file (300 DPI)
ggsave("ACME_ADE_Plot5.tiff", plot = acme_ade_plot5, dpi = 300, width = 8, height = 2, units = "in", device = "tiff", bg = "white")





###################################################################
###### Figure A9:  Causal mediation effects (the H&S measure only) 
###################################################################



## State capacity measure 1: the H&S measure only


# Mediator model (binary logistic regression)
mediator_model6 <- glm(bw_dummy ~ state_cap_HS + riv_dummy + threat + neighbors_number + prev_civ_conf + extek + mtnest + oilpcl + polity2_P4 + lnpop + postcold +
                         region_SSA + region_SoAsia + region_MENA + region_Latin + region_EAsiaPac + region_EurCentAsia + region_NoAmer, 
                       data = data, family = binomial)

summary(mediator_model6)



# Outcome model (binary logistic regression)
outcome_model6 <- glm(civ_war ~ state_cap_HS + bw_dummy + riv_dummy + threat + neighbors_number + prev_civ_conf + extek + +mtnest + oilpcl + polity2_P4 + lnpop + postcold +
                        region_SSA + region_SoAsia + region_MENA + region_Latin + region_EAsiaPac + region_EurCentAsia + region_NoAmer, 
                      data = data, family = binomial)

summary(outcome_model6)

# Perform mediation analysis
mediation_result <- mediate(mediator_model6, outcome_model6, treat = "state_cap_HS", mediator = "bw_dummy", boot = TRUE, sims = 1000)

# Summarize the mediation analysis results
summary(mediation_result)




### Figure 9 (Upper Panel)


# Tidy the mediator model coefficients
mediator_coeffs <- tidy(mediator_model6)

# Tidy the outcome model coefficients
outcome_coeffs <- tidy(outcome_model6)

# Combine the two datasets
mediator_coeffs$model <- "Mediator Model6"
outcome_coeffs$model <- "Outcome Model6"

coefficients <- bind_rows(mediator_coeffs, outcome_coeffs)

# Filter out the unwanted region variables and "mtnest"
coefficients_filtered <- subset(coefficients, 
                                !term %in% c("region_SSA", "region_SoAsia", "region_MENA", 
                                             "region_Latin", "region_EAsiaPac", 
                                             "region_EurCentAsia", "region_NoAmer", "mtnest"))

# Rename the terms, including renaming the intercept
coefficients_filtered$term <- recode(coefficients_filtered$term,
                                     "(Intercept)" = "Constant",
                                     "bw_dummy" = "Border barrier",
                                     "riv_dummy" = "Interstate rivalry",
                                     "threat" = "External threat",
                                     "neighbors_number" = "Number of neighbors",
                                     "prev_civ_conf" = "Number of previous civil wars",
                                     "extek" = "Excluded group with large TEK",
                                     "mtnest"= "Mountainous terrain",
                                     "oilpcl" = "Oil production per capita",
                                     "polity2_P4" = "Democracy",
                                     "lnpop" = "Total population (logged)",
                                     "state_cap_HS" = "The H&S Measure",
                                     "postcold" = "Cold War"
)

# Order the terms in the specified order (without Mountainous terrain)
order_of_terms <- c(
  "The H&S Measure", "Border barrier", "Interstate rivalry", "External threat", "Number of neighbors", 
  "Number of previous civil wars", "Excluded group with large TEK", "Mountainous terrain",
  "Oil production per capita", "Democracy", "Total population (logged)", 
  "Cold War", "Constant"
)
coefficients_filtered$term <- factor(coefficients_filtered$term, levels = rev(order_of_terms))

# Plot the filtered, renamed, and ordered coefficients with a reference line at x = 0
ggplot(coefficients_filtered, aes(x = estimate, y = term, color = model)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = estimate - std.error, xmax = estimate + std.error), height = 0.2) +
  facet_wrap(~model, scales = "free") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +  # Add reference line at x = 0
  labs(title = "",
       x = "Estimate",
       y = "Term") +
  theme_minimal()



# Save as high-resolution TIFF file (300 DPI, 2400x1200 pixels)
ggsave("Coefficient_Plot6.tiff", plot = last_plot(), dpi = 300, width = 8, height = 4, units = "in", device = "tiff", bg = "white")



### Figure 9 (Bottom Panel)


# Plot the ACME and ADE for the average effects with custom y-axis titles

acme_ade_average6 <- data.frame(
  effect = c("Average Causal Mediation Effect", "Average Direct Effect"),
  estimate = c(mediation_result$d.avg, mediation_result$z.avg),
  lower = c(mediation_result$d.avg.ci[1], mediation_result$z.avg.ci[1]),
  upper = c(mediation_result$d.avg.ci[2], mediation_result$z.avg.ci[2])
)


acme_ade_plot6 <- ggplot(acme_ade_average6, aes(x = estimate, y = effect)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +  # Add reference line at x = 0
  labs(title = "",
       x = "Estimate",
       y = "") +
  theme_minimal()+
  theme(
    panel.background = element_rect(fill = "white", color = NA), # Ensure white background
    plot.background = element_rect(fill = "white", color = NA)
  )

# Save as high-resolution TIFF file (300 DPI)
ggsave("ACME_ADE_Plot6.tiff", plot = acme_ade_plot6, dpi = 300, width = 8, height = 2, units = "in", device = "tiff", bg = "white")








